%the revisible part are K,N,RNK, R(s), sigma(s), y0, tspan
clear all
% parametor of index J
K=30; N=4; RNK=[4,4,4,3,3,3,2,2,2,2,ones(1,20)];%,1,1,1,1,1,1,1,1,1,1]; %control index
K=70; N=1; RNK=ones(1,70);
% K=3; N=4; RNK=[4,2,1];
% construct JNKR
JNKR = index(K,N,RNK);
fname='indexJNKR';
save(fname,'JNKR');
load(fname,'JNKR');
% time 
t = 1/2;
% construct orthogoal basis in L2[0,t]
orthbasis(t);
% interest
R=@(s) 0+s-s;
% volatility

% initial condition
y0=zeros(size(JNKR,1),1);
yic=0;
y0(1)=yic;
% solver
tspan = linspace(0,t,51);
options = odeset('AbsTol',1e-12,'RelTol',1e-6);
%for mean
%%Y=zeros(51,1);
%for var
Y=zeros(51,size(JNKR,1));

numpath=1000;
numpath=1;
for ii=1:numpath
    disp(ii)
    sig=ouprocess(0,0,0,1,5001);
    sig=2*ones(size(sig));
    %sig=exp(linspace(0,1/2,5001));
    %sig=sin(linspace(0,1/2,5001));
    [time,y] = ode45(@randomforce,tspan,y0,[],R,sig);
    %Y=Y+y(:,1); %for mean
    Y=Y+y.*y; %for var
end
% error
%meanerr(time',y,R,sig,y0,Y/numpath)
varerr(time',y,R,sig,y0,Y/numpath)


